﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using MathNet.Numerics.LinearAlgebra.Double;
using System.IO;
using System.Diagnostics;
using System.Threading.Tasks;
using System.Drawing;
using System.Windows.Forms;
using MathNet.Numerics.LinearAlgebra.Generic;
using System.Drawing.Imaging;

namespace T3Kmean
{
    /// <summary>
    /// 
    /// </summary>
    public interface IClusterMix
    {
        /// <summary>
        /// Updates this Model.
        /// </summary>
        void Update();
        /// <summary>
        /// Checks the errors of the points in this model.
        /// </summary>
        /// <returns></returns>
        Tuple<int, int, int>[] Check();

        /// <summary>
        /// Gets the R.
        /// </summary>
        double R { get; }
        /// <summary>
        /// Gets a value indicating whether this <see cref="IClusterMix"/> is modified in the last update.
        /// </summary>
        /// <value>
        ///   <c>true</c> if modified; otherwise, <c>false</c>.
        /// </value>
        bool Modified { get; }

    }

    public class Point
    {
        private Cluster cluster;
        private DenseVector bits;
        private int label;

        /// <summary>
        /// Gets or sets the cluster where the point is.
        /// </summary>
        /// <value>
        /// The cluster.
        /// </value>
        public Cluster Cluster
        {
            get { return cluster; }
            set {  }
        }
        /// <summary>
        /// Gets or sets the bits of the point.
        /// </summary>
        /// <value>
        /// The bits.
        /// </value>
        public DenseVector Bits
        {
            get { return bits; }
            private set { bits = value; }
        }
        /// <summary>
        /// Gets the label of the point.
        /// </summary>
        public int Label
        {
            get { return label; }
            private set { }
        }


        private double norm;

        /// <summary>
        /// Initializes a new instance of the <see cref="Point"/> class.
        /// </summary>
        /// <param name="s">The line of the file.</param>
        public Point(string s)
        {
            var data = s.Split(new char[] { ' ' });
            var _bits = (from b in data.Take(256)
                         select double.Parse(b)).ToArray();
            bits = new DenseVector(_bits);
            var lbl = data.Skip(256).ToArray();
            for (int i = 0; i < lbl.Length; i++)
            {
                if (lbl[i] == "1")
                    label = i;
            }
        }
        /// <summary>
        /// Sets the cluster where the point is.
        /// </summary>
        /// <param name="c">The cluster.</param>
        public void SetCluster(Cluster c)
        {
            c.AddPoint(this);
            norm = c.GetSQDistance(this);
            cluster = c;
        }
        /// <summary>
        /// Gets the gamma of the point.
        /// </summary>
        /// <param name="mu">The mu of the cluster.</param>
        /// <returns></returns>
        public double GetBernoulliProb(DenseVector mu)
        {
            double res = 1;
            for (int i = 0; i < mu.Count; i++)
            {
                res *= (bits[i] == 0 ? 1 - mu[i] : mu[i]);
            }
            return res;
        }
    }

    public class Cluster
    {
        public static DenseVector zero = new DenseVector(256);
        private DenseVector mu;
        private List<Point> points = new List<Point>();
        /// <summary>
        /// Gets the points in the cluster.
        /// </summary>
        public List<Point> Points
        {
            get { return points; }
            private set { }
        }
        /// <summary>
        /// Gets or sets the mu of the cluster.
        /// </summary>
        /// <value>
        /// The mu.
        /// </value>
        public DenseVector Mu
        {
            get { return mu; }
            set { mu = value; }
        }

        /// <summary>
        /// Create a new cluster froms the points.
        /// </summary>
        /// <param name="pts">A list of points.</param>
        /// <returns></returns>
        public static Cluster FromPoint(List<Point> pts)
        {
            Cluster c = new Cluster(Cluster.zero);
            for (int i = 0; i < pts.Count; i++)
            {
                pts[i].SetCluster(c);
            }
            c.UpdateCentroid();
            return c;
        }
        /// <summary>
        /// Initializes a new instance of the <see cref="Cluster"/> class.
        /// </summary>
        /// <param name="initial">The initial mu.</param>
        public Cluster(DenseVector initial)
        {
            mu = initial;
        }
        /// <summary>
        /// Adds the point.
        /// </summary>
        /// <param name="p">The point.</param>
        public void AddPoint(Point p)
        {
            points.Add(p);
        }
        /// <summary>
        /// Update the centroid of the cluster.
        /// </summary>
        public void UpdateCentroid()
        {
            if (points.Count != 0)
                mu = points.Aggregate<Point, DenseVector>(zero, (acc, p) => acc + p.Bits) / points.Count;
        }
        /// <summary>
        /// Clears the points.
        /// </summary>
        public void ResetPoints()
        {
            points.Clear();
        }
        /// <summary>
        /// Gets the Squared distance between the point and the centroid of the cluster.
        /// </summary>
        /// <param name="p">The point.</param>
        /// <returns></returns>
        public double GetSQDistance(Point p)
        {
            return (p.Bits - mu).Aggregate(0.0, (acc, d) => acc + (d * d));
        }
        /// <summary>
        /// Gets the dominants label of the points in the clusters.
        /// </summary>
        /// <returns></returns>
        public int DominantLabel()
        {
            var counts = new MathNet.Numerics.LinearAlgebra.Single.DenseVector(10);

            for (int i = 0; i < points.Count; i++)
            {
                counts[points[i].Label]++;
            }
            return counts.MaximumIndex();
        }
    }

    public class K_MeansModel : IClusterMix
    {
        private List<Point> points;
        private Cluster[] clusters;
        private bool modified = false;
        private Random seed = new Random();
        /// <summary>
        /// Gets a value indicating whether this <see cref="IClusterMix"/> is modified in the last update.
        /// </summary>
        /// <value>
        ///   <c>true</c> if modified; otherwise, <c>false</c>.
        /// </value>
        public bool Modified
        {
            get { return modified; }
            private set { }
        }
        /// <summary>
        /// Gets the clusters of the model.
        /// </summary>
        public Cluster[] Clusters
        {
            get { return clusters; }
            private set { }
        }

        //private K_MeansModel(int k, List<Point> pts)
        //{
        //    points = pts;
        //    clusters = new Cluster[k];
        //    var starting_mu = new DenseVector[k];
        //    for (int i = 0; i < k; i++)
        //    {
        //        var p = points[seed.Next(0, 1593)];
        //        if (starting_mu.Contains(p.Bits))
        //        {
        //            i--;
        //            continue;
        //        }
        //        else
        //        {
        //            starting_mu[i] = p.Bits;
        //        }
        //    }
        //    for (int i = 0; i < k; i++)
        //    {
        //        clusters[i] = new Cluster(starting_mu[i]);
        //    }
        //}
        /// <summary>
        /// Initializes a new instance of the <see cref="K_MeansModel"/> class.
        /// </summary>
        /// <param name="k">The number of clusters.</param>
        /// <param name="pts">The point list.</param>
        /// <param name="fake">if set to <c>true</c> [puts points with the same label in the same cluster].</param>
        public K_MeansModel(int k, List<Point> pts, bool fake)
        {
            if (fake)
            {
                points = pts;
                clusters = new Cluster[k];
                for (int i = 0; i < k; i++)
                {
                    var ps = pts.Where(p => p.Label == i).ToList();
                    clusters[i] = Cluster.FromPoint(ps);
                }
            }
            else
            {
                points = pts;
                clusters = new Cluster[k];
                var starting_mu = new DenseVector[k];
                for (int i = 0; i < k; i++)
                {
                    var p = points[seed.Next(0, 1593)];
                    if (starting_mu.Contains(p.Bits))
                    {
                        i--;
                        continue;
                    }
                    else
                    {
                        starting_mu[i] = p.Bits;
                    }
                }
                for (int i = 0; i < k; i++)
                {
                    clusters[i] = new Cluster(starting_mu[i]);
                }
            }
        }

        /// <summary>
        /// Assigns points in the most near cluster.
        /// </summary>
        private void AssignPoints()
        {
            modified = false;
            for (int i = 0; i < clusters.Length; i++)
            {
                clusters[i].ResetPoints();
            }

            for (int i = 0; i < points.Count; i++)
            {
                var idx = ArgMinCluster(points[i]);
                if (clusters[idx] != points[i].Cluster)
                {
                    modified = true;
                }
                points[i].SetCluster(clusters[idx]);
            }
        }
        /// <summary>
        /// Updates this Model.
        /// </summary>
        public void Update()
        {
            this.AssignPoints();
            if (modified)
            {
                Parallel.ForEach(clusters, c => c.UpdateCentroid());


                //foreach (var c in clusters)
                //{
                //    c.ReCentroid();
                //}
            }
        }
        /// <summary>
        /// Gets the min most near cluster.
        /// </summary>
        /// <param name="p">The point.</param>
        /// <returns></returns>
        private int ArgMinCluster(Point p)
        {
            int min_i = 0;
            double min_v = clusters[0].GetSQDistance(p);
            for (int i = 1; i < clusters.Length; i++)
            {
                var d = clusters[i].GetSQDistance(p);
                if (d < min_v)
                {
                    min_v = d;
                    min_i = i;
                }
            }
            return min_i;
        }
        /// <summary>
        /// Checks errors of the points in this model.
        /// </summary>
        /// <returns></returns>
        public Tuple<int, int, int>[] Check()
        {
            var errors = new Tuple<int, int, int>[clusters.Length];
            for (int c= 0; c < clusters.Length; c++)
			{
                var label = clusters[c].DominantLabel();
                int err = 0;
                foreach (var p in clusters[c].Points)
                {
                    if (p.Label != label)
                    {
                        err++;
                    }
                }
                errors[c] = new Tuple<int, int, int>(label, err, clusters[c].Points.Count);
                //Console.WriteLine("Cluster : " + label + " Errors : " + err);
            }
            return errors;
        }
        /// <summary>
        /// Gets the R.
        /// </summary>
        public double R
        {
            get
            {
                double a = 0, b = 0;
                for (int i = 0; i < points.Count; i++)
                {
                    for (int j = 0; j < points.Count; j++)
                    {
                        if (i != j)
                        {
                            var p1 = points[i];
                            var p2 = points[j];
                            if (p1.Label == p2.Label && p1.Cluster == p2.Cluster)
                            {
                                a++;
                            }
                            else if (p1.Label != p2.Label && p1.Cluster != p2.Cluster)
                            {
                                b++;
                            }
                        }
                    }
                }
                return 2 * (a + b) / ((double)(points.Count * (points.Count - 1)));
            }
        }
    }

    public class BernulliMix : IClusterMix
    {
        private ParallelOptions p_o = new ParallelOptions();
        private List<Tuple<Point, DenseVector>> points;
        private int k;
        private DenseVector pis;
        private DenseMatrix clusters;
        private bool modified = false;
        private Random seed = new Random();
        private double threshold;
        /// <summary>
        /// Gets a value indicating whether this <see cref="IClusterMix"/> is modified in the last update.
        /// </summary>
        /// <value>
        ///   <c>true</c> if modified; otherwise, <c>false</c>.
        /// </value>
        public bool Modified
        {
            get { return modified; }
            private set { }
        }
        /// <summary>
        /// Gets the matrix containing the mu of the clusters.
        /// </summary>
        public DenseMatrix Clusters
        {
            get { return clusters; }
            private set { }
        }

        /// <summary>
        /// Computes the gammas.
        /// </summary>
        private void ComputeGammas()
        {
            Parallel.For(0, points.Count, p_o, GammaBody);
            //for (int i = 0; i < points.Count; i++)
            //{
            //    GammaBody(i);
            //}
        }

        /// <summary>
        /// The body of the gamma for cycle.
        /// </summary>
        /// <param name="i">The i.</param>
        private void GammaBody(int i)
        {
            var p = points[i].Item1;
            var down = 0.0;
            var probs = new double[clusters.RowCount];
            for (int j = 0; j < clusters.RowCount; j++)
            {
                probs[j] = pis[j] * p.GetBernoulliProb(clusters.Row(j) as DenseVector);
                down += probs[j];
            }
            for (int j = 0; j < clusters.RowCount; j++)
            {
                points[i].Item2[j] = probs[j] / down;
            }
        }
        /// <summary>
        /// Computes the mu of clusters of the model.
        /// </summary>
        private void ComputeMu()
        {
            var ns = new DenseVector(clusters.RowCount);
            Parallel.For(0, clusters.RowCount, p_o, i => ns[i] = points.Sum(p => p.Item2[i]));
            //for (int i = 0; i < clusters.RowCount; i++)
            //{
            //    ns[i] = points.Sum(p => p.Item2[i]);
            //}
            

            for (int i = 0; i < clusters.RowCount; i++)
            {
                Vector<double> a = 1.0 / ns[i] * points.Aggregate(new DenseVector(points[0].Item1.Bits.Count), (acc, pg) => acc + (pg.Item2[i] * pg.Item1.Bits));
                
                if (!modified)
                {
                    for (int j = 0; j < a.Count; j++)
                    {
                        if (Math.Abs(a[j] - clusters.Row(i)[j]) >= threshold)
                        {
                            modified = true;
                            break;
                        }
                    }
                }
                clusters.SetRow(i, a);
            }
        }
        /// <summary>
        /// Computes the PI vector of the model.
        /// </summary>
        private void ComputePI()
        {
            pis = points.Aggregate(new DenseVector(k), (acc, p) => acc + p.Item2) / (double)points.Count;
        }


        /// <summary>
        /// Updates this Model.
        /// </summary>
        public void Update()
        {
            modified = false;
            ComputeGammas();
            ComputeMu();
            ComputePI();
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="BernulliMix"/> class.
        /// </summary>
        /// <param name="_points">The points.</param>
        /// <param name="_threshold">The threshold.</param>
        /// <param name="count">The number of cluster.</param>
        public BernulliMix(List<Point> _points, double _threshold, int count)
        {
            p_o.MaxDegreeOfParallelism = 8;
            threshold = _threshold;
            k = count;
            points = new List<Tuple<Point, DenseVector>>();
            foreach (var p in _points)
            {
                points.Add(new Tuple<Point, DenseVector>(p, new DenseVector(k)));
            }
            var ps = (from i in Enumerable.Range(0, k) select 1.0 / (double)k);
            pis = new DenseVector(ps.ToArray());
            clusters = new DenseMatrix(k, points[0].Item1.Bits.Count);
            for (int i = 0; i < clusters.RowCount; i++)
            {
                var mu = (from j in Enumerable.Range(0, points[0].Item1.Bits.Count) select (seed.NextDouble() * 0.5) + 0.25);
                clusters.SetRow(i, mu.ToArray());
            }
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="BernulliMix"/> class.
        /// </summary>
        /// <param name="_points">The _points.</param>
        /// <param name="_threshold">The _threshold.</param>
        /// <param name="count">The count.</param>
        public BernulliMix(K_MeansModel m, List<Point> _points, double _threshold)
        {
            p_o.MaxDegreeOfParallelism = 8;
            threshold = _threshold;
            k = m.Clusters.Length;
            points = new List<Tuple<Point, DenseVector>>();
            foreach (var p in _points)
            {
                points.Add(new Tuple<Point, DenseVector>(p, new DenseVector(k)));
            }
            var ps = (from i in Enumerable.Range(0, k) select 1.0 / (double)k);
            pis = new DenseVector(ps.ToArray());
            clusters = new DenseMatrix(k, points[0].Item1.Bits.Count);
            for (int i = 0; i < clusters.RowCount; i++)
            {
                Vector<double> mu = m.Clusters[i].Mu;
                clusters.SetRow(i, mu);
            }
        }

        /// <summary>
        /// Checks errors of the points in this model.
        /// </summary>
        /// <returns></returns>
        public Tuple<int, int, int>[] Check()
        {
            var label =
                (from p in points
                select new Tuple<Point, int>(p.Item1, p.Item2.MaximumIndex())).ToArray();

            var errors = new Tuple<int, int, int>[clusters.RowCount];
            var counts = (from i in Enumerable.Range(0, clusters.RowCount) select new DenseVector(10)).ToArray();


            for (int i = 0; i < label.Length; i++)
            {
                counts[label[i].Item2][points[i].Item1.Label]++;
            }


            for (int c = 0; c < clusters.RowCount; c++)
            {
                var dominant = counts[c].MaximumIndex();
                var inside = 0;

                int err = 0;
                foreach (var p in label)
                {
                    if (p.Item2 == c)
                    {
                        inside++;
                        if (p.Item1.Label != dominant)
                        {
                            err++;
                        }
                    }
                }
                errors[c] = new Tuple<int, int, int>(dominant, err,inside);
                //Console.WriteLine("Cluster : " + label + " Errors : " + err);
            }
            return errors;


        }

        /// <summary>
        /// Gets the R.
        /// </summary>
        public double R
        {
            get
            {
                double a = 0, b = 0;
                for (int i = 0; i < points.Count; i++)
                {
                    for (int j = 0; j < points.Count; j++)
                    {
                        if (i != j)
                        {
                            var p1 = points[i].Item1;
                            var c1 = points[i].Item2.MaximumIndex();
                            var p2 = points[j].Item1;
                            var c2 = points[j].Item2.MaximumIndex();
                            if (p1.Label == p2.Label && c1 == c2)
                            {
                                a++;
                            }
                            else if (p1.Label != p2.Label && c1 != c2)
                            {
                                b++;
                            }
                        }
                    }
                }
                return 2 * (a + b) / ((double)(points.Count * (points.Count - 1)));
            }
        }

        //gamma, media, pi;
        //medie tra 0.25 e 0.75
        //pi = 1/k
    }

    public class K_means_bernoulli : IClusterMix
    {
        private enum ModelType
        {
            BernoulliMixture,
            KMeans
        }
        public IClusterMix m;
        ModelType actual;
        private bool modified = false;
        private List<Point> points;
        private double threshold;

        /// <summary>
        /// Initializes a new instance of the <see cref="K_means_bernoulli"/> class.
        /// </summary>
        /// <param name="_points">The points.</param>
        /// <param name="_threshold">The threshold.</param>
        /// <param name="k">The k.</param>
        public K_means_bernoulli(List<Point> _points, double _threshold,  int k)
        {
            m = new K_MeansModel(k, _points, false);
            points = _points;
            actual = ModelType.KMeans;
            threshold = _threshold;
        }

        /// <summary>
        /// Updates this Model.
        /// </summary>
        public void Update()
        {
            switch (actual)
            {
                case ModelType.BernoulliMixture:
                    m.Update();
                    modified = m.Modified;
                    break;
                case ModelType.KMeans:
                    m.Update();
                    if (!m.Modified)
                    {
                        m = new BernulliMix(m as K_MeansModel, points, threshold);
                        actual = ModelType.BernoulliMixture;
                        Console.WriteLine("changedchangedchangedchangedchangedchangedchangedchangedchangedchanged");
                    }
                    modified = true;
                    break;
                default:
                    break;
            }
        }

        /// <summary>
        /// Checks the errors of the points in this model.
        /// </summary>
        /// <returns></returns>
        public Tuple<int, int, int>[] Check()
        {
            return new Tuple<int, int, int>[1];
        }

        /// <summary>
        /// Gets the R.
        /// </summary>
        public double R
        {
            get { return m.R; }
        }

        /// <summary>
        /// Gets a value indicating whether this <see cref="IClusterMix"/> is modified in the last update.
        /// </summary>
        /// <value>
        ///   <c>true</c> if modified; otherwise, <c>false</c>.
        /// </value>
        public bool Modified
        {
            get { return modified; }
        }
    }

    static class Program
    {
        static Stopwatch chrono = new Stopwatch();
        static Graphics graphics;

        private static DenseMatrix ToSQMatrix(this DenseVector v)
        {
            var res = new DenseMatrix(16);
            for (int i = 0; i < v.Count; i++)
            {
                res[i % 16, i / 16] = v[i];
            }

            return res;
        }

        private static void DrawClusters(Cluster[] cluster, DenseMatrix[] old)
        {
            for (int idx = 0; idx < cluster.Length; idx++)
            {
                DenseMatrix m = cluster[idx].Mu.ToSQMatrix();
                for (int j = 0; j < m.RowCount; j++)
                {
                    for (int i = 0; i < m.ColumnCount; i++)
                    {
                        var r = (m[i, j] != old[idx][i, j]) ? (int)(m[i, j] * 255.0) : 0;
                        var b = (m[i, j] == old[idx][i, j]) ? (int)(m[i, j] * 255.0) : 0;
                        var brush = new SolidBrush(Color.FromArgb(255, r, 0, b));
                        graphics.FillRectangle(brush, new Rectangle((i * 8) + (idx % 4) * 150, (j * 8) + 150 * (idx / 4), 8, 8));
                    }
                }
            }
        }

        private static void DrawBernoulli(DenseMatrix cluster, DenseMatrix[] old)
        {
            for (int idx = 0; idx < cluster.RowCount; idx++)
            {
                DenseMatrix m = (cluster.Row(idx) as DenseVector).ToSQMatrix();
                for (int j = 0; j < m.RowCount; j++)
                {
                    for (int i = 0; i < m.ColumnCount; i++)
                    {
                        var cond = Math.Abs(m[i, j] - old[idx][i, j]) < 0.01;
                        var r = (!cond) ? (int)(m[i, j] * 255.0) : 0;
                        var b = (cond) ? (int)(m[i, j] * 255.0) : 0;
                        var brush = new SolidBrush(Color.FromArgb(255, r, 0, b));
                        graphics.FillRectangle(brush, new Rectangle((i * 8) + (idx % 4) * 150, (j * 8) + 150 * (idx / 4), 8, 8));
                        
                    }
                }
            }
        }

        private static void SaveBernoulli(DenseMatrix cluster, Graphics gr, int voffset) //, string path)
        {
            //var temp = new Bitmap(49 * cluster.RowCount, 48);
            //var gr = Graphics.FromImage(temp);

            for (int idx = 0; idx < cluster.RowCount; idx++)
            {
                DenseMatrix m = (cluster.Row(idx) as DenseVector).ToSQMatrix();
                for (int j = 0; j < m.RowCount; j++)
                {
                    for (int i = 0; i < m.ColumnCount; i++)
                    {
                        var brush = new SolidBrush(Color.FromArgb(255, 0, 0, (int)(m[i, j] * 255.0)));
                        gr.FillRectangle(brush, new Rectangle((i * 3) + idx * 49, (j * 3) + (voffset * 49), 3, 3));

                    }
                }
            }


            //temp.Save(path, ImageFormat.Png);
        }

        private static void SaveKMeans(Cluster[] cluster, Graphics gr, int voffset) //, string path)
        {
            //var temp = new Bitmap(49 * cluster.Length, 48);
            //var gr = Graphics.FromImage(temp);

            for (int idx = 0; idx < cluster.Length; idx++)
            {
                DenseMatrix m = cluster[idx].Mu.ToSQMatrix();
                for (int j = 0; j < m.RowCount; j++)
                {
                    for (int i = 0; i < m.ColumnCount; i++)
                    {
                        var brush = new SolidBrush(Color.FromArgb(255, 0, 0, (int)(m[i, j] * 255.0)));
                        gr.FillRectangle(brush, new Rectangle((i * 3) + idx * 49, (j * 3) + (voffset * 49), 3, 3));
                    }
                }
            }


            //temp.Save(path, ImageFormat.Png);
        }

        public static void SaveAll(K_MeansModel m1, BernulliMix m2, K_means_bernoulli m12, string path)
        {
            var temp = new Bitmap(49 * m1.Clusters.Length, (49 * 3));
            var gr = Graphics.FromImage(temp);

            SaveKMeans(m1.Clusters, gr,0);
            SaveBernoulli(m2.Clusters, gr, 1);
            SaveBernoulli((m12.m as BernulliMix).Clusters , gr,2);

            temp.Save(path, ImageFormat.Png);
        }

        private enum ClusterType
        {
            BernoulliMixture,
            KMeans
        }

        static void Main(string[] args)
        {
            System.Threading.Thread.CurrentThread.CurrentCulture = System.Globalization.CultureInfo.CreateSpecificCulture("en");
            System.Threading.Thread.CurrentThread.CurrentUICulture = System.Globalization.CultureInfo.CreateSpecificCulture("en");

            var str = File.CreateText("out.txt");
            
            var text = File.ReadAllLines("semeion.data");
            var points =
                   (
                       from line in text
                       where line != ""
                       select new Point(line)
                   ).ToList();

            Console.WriteLine("Show images? (yes/no)");
            bool show = Console.ReadLine() != "no";

            int k = -1;
            Console.WriteLine("How many Clusters? (min = 5; default = 10; max = 15)");
            if (int.TryParse(Console.ReadLine(), out k))
                k = Clamp(k, 5, 15);
            else 
                k = 10;

            Console.WriteLine("Fake initialization? (yes/no)");
            bool fake = Console.ReadLine() == "yes";

            Console.WriteLine("K-Means or Bernoulli mixture? (1/2)");
            ClusterType type = Console.ReadLine() == "1" ? ClusterType.KMeans : ClusterType.BernoulliMixture;
            IClusterMix m = new K_MeansModel(0, new List<Point>(), false);
            switch (type){
                case ClusterType.BernoulliMixture:
                    m = new BernulliMix(points, 0.01, k);
                    break;
                case ClusterType.KMeans:
                    m = new K_MeansModel(k, points, fake);
                    break;
            }
            Form f = null;
            if (show)
            {
                f = new Form();
                f.Width = 600;
                f.Height = 650;
                f.StartPosition = FormStartPosition.CenterScreen;
                graphics = f.CreateGraphics();
                f.Show();
            }

            for (int ks = 5; ks <= 15; ks++)
            {                
                
                var iterations = new double[10];
                var rs = new double[10];
                var errors = new List<Tuple<int, int, int>>[10];
                var times = new List<double>();
                var drs = new double[10];
                for (int i = 0; i < 1; i++)
                {
                    var m1 = new K_MeansModel(ks, points, false);
                    var m2 = new BernulliMix(points, 0.01, ks);
                    var m12 = new K_means_bernoulli(points, 0.01, ks);

                    //m = new K_means_bernoulli(points, 0.01, ks);
                    ////m = new K_MeansModel(ks, points, false);
                    //int iteration = 0;
                    //do
                    //{
                    //    DenseMatrix[] old = new DenseMatrix[1];
                    //    if (show)
                    //    {
                    //        switch (type)
                    //        {
                    //            case ClusterType.KMeans:
                    //                old =
                    //                   (from c in (m as K_MeansModel).Clusters
                    //                    select c.Mu.ToSQMatrix()).ToArray();
                    //                break;
                    //            case ClusterType.BernoulliMixture:
                    //                old =
                    //               (from c in Enumerable.Range(0, (m as BernulliMix).Clusters.RowCount)
                    //                select ((m as BernulliMix).Clusters.Row(c) as DenseVector).ToSQMatrix()).ToArray();
                    //                break;
                    //        }
                    //    }

                    //    //var old_r = m.R;
                    //    chrono.Restart();
                    //    m.Update();
                    //    chrono.Stop();
                    //    //var r = m.R;
                    //    Console.WriteLine(iteration++);
                    //    //Console.WriteLine(iteration++ + " in " + chrono.Elapsed.TotalMilliseconds.ToString("0.00") + " --> " + r + " --> " + Math.Sign(r - old_r));
                    //    iterations[i]++;
                    //    times.Add(chrono.Elapsed.TotalMilliseconds);
                    //    //rs[i] = r;
                    //    //drs[i] = Math.Sign(r - old_r);

                    //    if (show)
                    //    {
                    //        graphics.Clear(Color.White);
                    //        switch (type)
                    //        {
                    //            case ClusterType.KMeans:
                    //                DrawClusters((m as K_MeansModel).Clusters, old);
                    //                break;
                    //            case ClusterType.BernoulliMixture:
                    //                DrawBernoulli((m as BernulliMix).Clusters, old);
                    //                break;
                    //        }
                    //        f.Update();
                    //    }
                    //} while (m.Modified);

                    do
                    {
                        m1.Update();
                    } while (m1.Modified);

                    do
                    {
                        m2.Update();
                    } while (m2.Modified);

                    do
                    {
                        m12.Update();
                    } while (m12.Modified);

                    var act_errors = m.Check();
                    errors[i] = act_errors.ToList();
                    for (int j = 0; j < errors.Length; j++)
                    {
                        //Console.WriteLine("Cluster : " + j + " with dominant label : " + errors[j].Item1 + " Points : " + errors[j].Item3 + " errors : " + errors[j].Item2 + " rel error : " + ((double)errors[j].Item2 / (double)errors[j].Item3).ToString("0.00"));
                        //switch (type)
                        //{
                        //    case ClusterType.KMeans:
                        //        var km = m as K_MeansModel;
                        //        Console.WriteLine("Cluster : " + i + " Points : " + km.Clusters[i].Points.Count + " with dominant label : " + errors[i].Item1 + " errors : " + errors[i].Item2 + " rel error : " + ((double)errors[i].Item2 / (double)km.Clusters[i].Points.Count).ToString("0.00"));
                        //        break;
                        //    case ClusterType.BernoulliMixture:
                        //        Console.WriteLine("Cluster : " + i + " with dominant label : " + errors[i].Item1 + " errors : " + errors[i].Item2);
                        //        break;
                        //}
                    }
                    rs[i] = m.R;
                    Console.WriteLine(rs[i]);

                    SaveAll(m1, m2, m12, "all.png" + ks + ".png");
                    //SaveBernoulli((m as BernulliMix).Clusters, "Bernoulli" + ks + ".png");

                }
                    str.WriteLine("Bernoulli");
                    var ite = iterations.Mean();
                    var erre = rs.Mean();
                    var time = times.Mean();
                    var dr = drs.Mean();
                    str.WriteLine(ite + ";" + time.ToString("0.00") + ";" + erre + ";" + dr);
                    //str.WriteLine("\n\n errors");

                    
                    Console.WriteLine("ite ---- "+ks);
                    //for (int iter = 0; iter < k; iter++)
                    //{
                    //    var line = 
                    //        (from el in errors

                    //}

                
            }
            Console.WriteLine("R : " + m.R);
//            Console.WriteLine("R : " + iterations.Sum() / 10.0);
            Console.WriteLine("Premere un tasto per continuare . . . ");
            //Console.ReadKey();
            str.Close();
        }
        private static int Clamp(int i, int min, int max)
        {
            if (min <= i && max >= i)
                return i;
            if (i < min)
                return min;
            return max;
        }

        private static double Mean(this IEnumerable<double> x)
        {
            return x.Sum() / (double)x.Count();
        } 
    }
}
